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Abstract 

This work is devoted to the thermodynamics of high-temperature dense hy- 
drogen plasmas in the pressure region between 10~^ and 10^ Mbar. In par- 
ticular we present for this region results of extensive calculations based on 
a recently developed path integral Monte Carlo scheme (direct PIMC). This 
method allows for a correct treatment of the thermodynamic properties of 
hot dense Coulomb systems. Calculations were performed in a broad region 
of the nonideality parameter F < 3 and degeneracy parameter rigA^ < 10. 
We give a comparison with a few available results from other path integral 
calculations (restricted PIMC) and with analytical calculations based on Fade 
approximations for strongly ionized plasmas. Good agreement between the 
results obtained from the three independent methods is found. 
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I. INTRODUCTION 



The thermodynamics of strongly correlated Fermi systems at high pressures are of grow- 
ing importance in many fields, including shock and laser plasmas, astrophysics, solids and 
nuclear matter, see Refs. [0-0] for an overview. In particular the thermodynamical properties 
of hot dense plasmas under high pressure are of importace for the description of plasmas 
relevant for laser fusion Further among the phenomena of current interest are Fermi 
liquids, metallic hydrogen 0, plasma phase transition, e.g. ^ and references therein, bound 
states etc., which occur in situations where both Coulomb and quantum effects are relevant. 
There has been significant progress in recent years to study these systems analytically and 
numerically, see e.g. |'^- |T0| , p!3H l5| . Due to the enormeous difficulties to develop analytical 
descriptions for hydrogen plasmas with strong coupling, e.g. there is still an urgent 

need to test the analytical theory by an independent numerical approach. 

An approach which is particularly well suited to describe thermodynamic properties 
in the region of high pressure, characterized by strong coupling and strong degeneracy, 
is the path integral quantum Monte Carlo (PIMC) method. There has been remarkable 
recent progress in applying these techniques to Fermi systems, for an overview see e.g. 
Refs. [p!|p|,|T6|-[T8[] . However, these simulations are essentially hampered by the fermion sign 



problem. To overcome this difficulty, several strategies have been developed to simulate 
macroscopic Coulomb systems ||8| JT9|j20[| : the first is the restricted PIMC concept where 
additional assumptions on the density operator p are introduced which reduce the sum 
over permutations to even (positive) contributions only. This requires knowledge of the 
nodes of the density matrix which is available only in a few special cases, e.g. |ll9| , ^ . 
However, for interacting macroscopic systems, these nodes are known only approximately. 



e.g. |2T[, and the accuracy of the results is difficult to assess from within this scheme. An 
alternative are direct fermionic PIMC simulations which have occasionally been attempted 
by various groups but which were not sufficiently precise and efficient for practical 
purposes. Recently, three of us have proposed a new path integral representation for the N- 
particle density operator |23,24| which allows for direct fermionic path integral Monte Carlo 
simulations of dense plasmas in a broad range of densities and temperatures. Using this 
concept we computed the pressure and energy of a degenerate strongly coupled hydrogen 
plasma P^ , pB| and the pair distribution functions in the region of partial ionization and 
dissociation p6| , p7| . This scheme is rather efficient when the number of time slices (beads) 



in the path integral is less or equal 50 and was found to work well for temperatures ksT > 
O.lRy. In this paper we derive further improved formulas for the pressure and energy and 
give, for the first time, a detailed derivation of all main results and rigorously justify the 
use of the effective quantum pair potential (Kelbg potential) in direct PIMC simulations. 
Further, in the present work this method will be applied to high-pressure plasmas {p ~ 
10~^ — lO^Mbar) in such temperature regions were considerable deviations from the ideal 
behavior are observed. 

One difficulty of PIMC simulations is that reliable error estimates are often not available, 
in particular for strongly coupled degenerate systems. Moreover, in this region no reliable 
data from other theories such as density functional theory or quantum statistics, e.g. []3|,p!5 



are available which would allow for an unambiguous test. Furthermore, results from classi- 
cal molecular dynamics simulations exist, but they apply only to fully ionized and weakly 
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degenerate plasmas, e.g. p8|-|30|, which is outside the range of interest for this work. Also, 
new quantum molecular dynamics approaches are being developed, e.g. [|T0|-[l2[|, but they 
are only beginning to produce accurate results. 

Therefore, it is of high interest to perform quantitative comparisons of independent sim- 
ulations, such as restricted and direct fermionic PIMC, and to develop improved analytical 
approximations, which is the aim of this paper. We compare recent results of Militzer et al. 
3^ for pressure and energy isochors {n ~ 2.5- lO^^cm"^) of dense hydrogen to our own direct 



PIMC results. This is a non-trivial comparison since the two approaches employ indepen- 
dent sets of approximations. Nevertheless, we find very good agreement for temperatures 
ranging from lO^i^' to as low as 50, OOOi^. This is remarkable since there the coupling and 
degeneracy parameters reach rather large values, F ^ 3 and UgA^ ^ 10, and the plasma 
contains a substantial fraction of bound states. 

Further, we use the new data to make a comparison with analytical estimates which 
are based on Fade approximations for strongly ionized plasmas. These formulae were con- 
structed on the basis of the known analytical results for the limiting cases of low density 
PJ33|] and high density ||^. These Fade approximations are exact up to quadratic terms in 
the density and interpolate between the virial expansions and the high-density asymptotics 
P^^B[. We find that the results for the internal energy and for the pressure agree well 
with the FIMC results in the region of the density temperature plane, where F < 1.6 and 
nA^ < 5. 



II. PATH INTEGRAL REPRESENTATION OF THERMODYNAMIC 

QUANTITIES 

We now our direct FIMC scheme. All thermodynamic properties of a two-component 
plasma are defined by the partition function Z which, for the case of Ne electrons and Np 
protons, is given by 

with Q{Ne,Np,(3) = J2 dqdrp{q,r,a;l3), (1) 

- V 

where /3 = 1/kBT. The exact density matrix is, for a quantum system, in general, not known 



but can be constructed using a path integral representation [^| 



V " V 

X 



5^5^(±l)''-5(a,PaOPp(-+^), (2) 



a P 



where p» = p A/5) = (P^*"^) |e-^^^|i?»), whereas A/3 = [5/{n + 1) and AA^ = 

27r/i^A/3/ma, a = p,e. H is the Hamilton operator, H = K + tic-, containing kinetic and 
potential energy contributions, K and t/c, respectively, with Uc = + + U^^ being the 
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sum of the Coulomb potentials between protons (p), electrons (e) and electrons and protons 
(ep)]. Further, i?» = (gW,r«) = {R^J\r^:^), for z = 1, . . . n + 1, = (g,r) = 
and = and a' = a. This means, the particles are represented by fermionic 

loops with the coordinates (beads) [R] = [R^^^; R^^^; . . . ; i?'-"-*; R^^^^^], where q and r denote 
the electron and proton coordinates, respectively. The spin gives rise to the spin part 
of the density matrix S, whereas exchange effects are accounted for by the permutation 
operator P, which acts on the electron coordinates and spin projections, and the sum over 
the permutations with parity Kp. In the fermionic case (minus sign), the sum contains 
NJ/2 positive and negative terms leading to the notorious sign problem. Due to the large 
mass difference of electrons and ions, the exchange of the latter is not included. The matrix 
elements p*^*-* can be rewritten identically as 



To compute thermodynamic functions, the logarithm of the partition function has to be 
differentiated with respect to thermodynamic variables. In particular, for the equation of 
state p and internal energy E follows, 

(3p = d\nQ/dV = [a/3Vd\nQ/da]^=u (4) 
(3E = -(3d\nQ/d(3, (5) 

where a is a length scaling parameter a = L/Lq. This means, in the path integral repre- 
sentation (^, each high-temperature density matrix has to be differentiated in turn. For 
example, the result for the energy will have the form 



V 

n+l 

X 

k=l 



J]p«...p('=-i). .p('^+i)...p(")5^5^(±l)-5(a,Pa')Pp("+^), (6) 



cr P 



and, analogously for other thermodynamic functions. 

There are two different approaches to evaluate this expression. One is to first choose 
an approximation for the high-temperature density matrices p*^*^ and then to perform the 
differentiation. The other way is to first differentiate the operator expression for p^'^^ and 
use an approximation for the matrix elements only in the final result. As we checked, the 
second method is more accurate and will be used in the following. 

To evaluate the derivatives in Eq. @, it is convenient to indroduce dimensionless in- 
tegration variables t]'-'^-* = {7}^\r]^'')^ where 77!'^'' = K,a{R'a^ — R'a ioi k = 1, . . . ,n and 
a = p,e, and = maksT / {2Tifi?) = This has the advantage that now the dif- 



ferentiation of the density matrix affects only the interaction terms. Indeed, one can show 
that 
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where = {k,R^^\ k^R^^^), = (xf ,xi')), with xi'^ = KaR^^^ + Y.tiVa\ and k 
runs from 1 to n. Further, X*^"^-'^^ = {KpR^^^^\ Kg-Re"^"^'') = X*^"-', and we denoted 



g 47r(?i4-l) 
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e 2 



...ix^'^)), 



(8) 



where pi'^'' = pi'^V (i^aK), p{k) = {p^p\p'e^) and use has been made of Eq. (|]). For A; = n + 1, 
we have 



o- P 



cr P 



(9) 



Further, ?7c(X('=-^)) = U^c\x^^-^'^) + U^^\x^^-^^) , with f/i^^ and ^/P denoting the interac- 
tion between identical and different particle species, respectively, U^P{X) = U^{X) + Up{X) 
and UP{X) = U^P{X). 

Using these results and Eq. (0), we obtain for the energy 



3 

2' 



1 1 



i3E=-{N, + Npj ^ 3^ 



CT P 



X 



Epw.-.p^'^-^^ 

fc=i 

X p^'^) . . . p(")Pp("+^) 



P Tin P c\r, ' PPp 



dp 



dp 



X("+i)=X(0).ct'=ct 



(10) 



This way, the derivative of the density matrix has been calculated, and we turn to the next 
point - to find approximations for the high-temperature density matrix. 

III. HIGH-TEMPERATURE ASYMPTOTICS OF THE DENSITY MATRIX IN 
THE PATH INTEGRAL APPROACH. KELBG POTENTIAL 



In this section we derive an approximation for the high-tempature density matrix which 
is suitable for direct PIMC simulations. Further, we demonstrate that the proper choice 
of the effective quantum pair potential is given by the Kelbg potential. Following Refs. 
T6| , |38|j39| , we derive a modified representation for the density matrix. The mains steps are: 



1. The N-particle density matrix is expanded in terms of 2-particle, 3-particle etc. con- 
tributions from which only the first, Pab, is retained |lT6|,p8|,|39| ; 



2. In the high-temperature limit, pab factorizes into a kinetic (po) and an interaction term 
(Pu): Pab ~ PoPu 1 because it can be shown that 



I + O 



(11) 
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where / is the unity operator. In this way we get the following representation for the 
two-particle density matrix 

{mambf/^\ . nia , ,.21 r "^b / /x2l 



where $ab(ra, r'^, r;,, r'^) is the off-diagonal two-particle effective potential. 

3. In the following, the off-diagonal matrix elements of the effective binary potentials 
will be approximated by the diagonal ones by taking the Kelbg potential at the center 
coordinate, <l>'*^(r, r'; A/5) ^ A/3).; 

4. For the plasma parameter region of interest, the protons can be treated classically, and 

may be approximated by the Coulomb potential. 

We will now comment on these steps in some more detail. We calculated the effective 
potential by solving a Bloch equation by first order perturbation theory. The method has 
been described in detail in This procedure defines an effective off-diagonal quantum 
pair potential for Coulomb systems, which depends on the inter-particle distances Tab, r'ab- 
As a result of first-order perturbation theory we get explicitely 

$^^(r.„ rl„ A/3) ^ e.e, /' -^erf ( —^j^^=^] , (13) 

Jo dab[a) \2\ab^ a{l - a) ) 

where dab{a) = \avab + (1 — '^)''a6l' erf(x) is the error function erf(a;) = dte~^^ , and 

^ab — \^ with = +m^^. It is interesting to note, that a simple approximation of 
the complicated integral over a by the length of the interval multiplied with the integrand in 
the center (Mittelwertsatz) leads us to the so-called KTR-potential due to Klakow, Toepffer 
and Reinhard which (in the diagonal approximation) is often used in quasi-classical MD 



simulations 10,13 



In our direct PIMC calculations we used the full expression for the interaction potential, 
keeping the a-integration but, in order to save computer time, we approximated the two- 
particle interaction potential by its diagonal elements. The diagonal element (r'^^ = Vab) 
of is just the familiar Kelbg potential, given by (we will use the same notation for the 
potential) 

^'^^(Ir,,!, A/?) = <l>'^^(r,,, r,,. A/?) = [l - e'^^ + v^7r}x„5 (1 - erf(x„5))l ' ^^^^ 

'^ab-^ab '- 

where Xab = l^abl/^ab, and we underline that the Kelbg potential is finite at zero distance. 
The error of the above approximations, for each of the high-temperature factors on the r.h.s. 
of Eq. (H), is of the order l/{n + 1)^ 



With these approximations, we obtain the result p^'^^ = p'^^ pl'l^ + 0[{l/n + iy], where 



Po 



(k) 



is the kinetic density matrix, while p^^ = e ^P^i^''^ ^^^6{X'^^ — X'^^^), where U denotes 
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the followine; sum of Coulomb and Kelbg potentials, f/(X(^)) = U^{X^^^] 
U^P{X^''\x^e''^). Notice that special care has to be taken in performing the derivatives 
with respect to /? of the Coulomb potentials which appear in Eq. (|iy). Indeed, products 

p^-^^ . . . _Pp("+i) p ^^^'^"q'^ have a singularity at zero interparticle distance which is 

integrable but leads to difficulties in the simulations. Due to the Kelbg potential, for the e-e 
and p-p interaction, this singularity is weakenend, but it is enhanced for the e-p interaction. 
In order to assure efficient simulations we, therefore, further transform the e-p contribution 
in the following way: 



da 



g-A/3(l-a)i-|^(fc)^ 



(A/3t/P(i?(^-^))) 



-Af3K 



d 



-(3^ (A/5f/(2)(/?('=-i))) 



+ 0[(l/n + l)2]. 



(16) 



This means, within the standard error of our approximation O [(A/?)^], we have replaced 
the e-p Coulomb potential f/i^"* by the corresponding Kelbg potential f/*^^-*, which is much 
better suited for MC simulations. 

Thus, using Ap ^ Ag, we finally obtain for the energy: 



1 



1 



Q Ar^AAP^ 



r A^, 



dqdrd^Psiq, [r],/5) x 

AT 



1=0 I- p<t 



+E 



\' pt\ 



+ EE*? 

p=l t=l 



ep 



p<t 



+EE« 



p=i t=i - I P* I 
1 5det|C.^s 



n,l| 



detlCb^l 



with 



'^pt 



KtWpt) 



^pt 



dp 



2\xi 



pt\ 



(17) 



and = ApdlP'^^PilxLl, l3')]/dp'\i3'=A(3 contains the electron-proton Kelbg potential ^"p. 



Here, (. . . | . . .) denotes the scalar product, and g^f, rpt and Xpt are differences of two coor- 
dinate vectors: qpt = qp - qt, Vpt = Vp - rt, Xpt = Vp - qt, r^^ = Vpt + y^^, Xp^ = Xpt + Up and 

Upt = yp~ Vti with = AAg X]fc=i ^a'^ ■ Here we introduced dimensionless distances between 
neighboring vertices on the loop, C,''^\ . . . thus, explicitly, [r] = [r; yi^^; yi^^; . . .]. Further, 
the density matrix ps in Eq. (^) is given by 



n Ne 



Ne 



-(3U{q,[r],l3) 



71,1 I 

ab ls) 



(18) 



1=1 
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where U{q, [r], (3) = f/P(g) + [r], A/5) + f/^P(g, [r], A/3)}/(n + 1) and (^^^ ^ exp[-7r|e?^|2] 



We underhne that the density matrix (|TBp does not contain an exphcit sum over the per- 
mutations and thus no sum of terms with alternating sign. Instead, the whole exchange 
problem is contained in a single exchange matrix given by 



n,l| 
ab I 



r\{ra-ri,)+y^f 



(19) 



As a result of the spin summation, the matrix carries a subscript s denoting the number of 



electrons having the same spin projection. For more details, we refer to Refs. |[23|j24 
In similar way, we obtain the result for the equation of state. 



ppV 



1 + 



1 



Ne + Np A^'AAP' 



e s=Q 



X 



EBe^ Ape 



I dqdrdips{q, [r],/?) 

dA(3^^p 



pt 



2 

p=i t=i 



1=1 



d\xpt\ 
dAf^^^p 



.p<t 



pt\ 



p=i t=i 



dWpt\ 



+ 



with A, 



pt 



Kt 


rpt) 




"^pt 





a 




det|Cb^|« 


da 


^pt 


_ {^pt\^pt) 
\-^pt\ 



(20) 



The structure of Eqs. (|T7|, ^D]) is obvious: we have separated the classical ideal gas part 
(first term). The ideal quantum part in excess of the classical one and the correlation 
contributions are contained in the integral term, where the second line results from the ionic 
correlations (first term) and the e-e and e-i interaction at the first vertex (second and third 
terms respectively). The third and fourth lines are due to the further electronic vertices and 
the explicit temperature dependence [in Eq. (|T^ and volume dependence in Eq. (PP])] of the 
exchange matrix, respectively. The main advantage of Eqs. ([T7| , pOf ) is that the explicit sum 
over permutations has been converted into the spin determinant which can be computed 
very efficiently using standard linear algebra methods. Furthermore, each of the sums in 
curly brackets in Eqs. ([l^, ^) is bounded as the number of vertices increases, n — oo, and is 
thus well suited for efficient Monte Carlo simulations. Notice also that Eqs. ([l^, ^0|) contain 
the important limit of an ideal quantum plasma in a natural way [H^ . 



IV. COMPARISON OF DIRECT AND RESTRICTED PIMC SIMULATIONS 

Expressions (|1^ ^0]) are well suited for numerical evaluation using Monte Carlo tech- 
niques, e.g. |]T^|T^. In our Monte Carlo scheme we used three types of steps, where either 
electron or proton coordinates, or gj or inidividual electronic beads ^l'^^ were moved un- 
til convergence of the calculated values was reached. Our procedure has been extensively 
tested. In particular, we found from comparison with the known analytical expressions for 
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pressure and energy of an ideal Fermi gas that the Fermi statistics is very well reproduced 
26| . Further, we performed extensive tests for few-electron systems in a harmonic trap 



where, again, the analytically known limiting behavior (e.g. energies) is well reproduced 
43| , |4^ . For the present simulations of dense hydrogen, we varied both the particle number 



and the number of time slices (beads). As a result of these tests, we found that to obtain 
convergent results for the thermodynamic properties of dense hydrogen, particle numbers 
Ne = Np = 50 and beads numbers in the range of n = 6 ... 20 are adequate ||2 
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FIGURES 




log T [K] 



FIG. 1. Density-temperature plane showing the parameter region for which calculations are 
performed. The data of Fig. 2 are along the dashed line (isochor = 1.86). The data of Figs. 3 
and 4 are inside the bold rhomb, along lines of constant F between the lines nA^ = 2 and nA^ = 5, 
respectively. Data for the vertical line (isotherm T = 50, OOOK) are given in Fig. 5. 

We will now compare our results with some available results obtained by the Monte 



Carlo technique developed by the Urbana group ||19| , |32| . We may first state that both 
Monte Carlo techniques differ in several fundamental points, so that they are essentially 
independent approaches. Let us briefly outline the main differences between the technique 
developed in Urbana, known as the restricted PIMC scheme and references therein, and 
the approach described here. These authors performed simulations with 32 electrons and 
protons; their restricted PIMC scheme required to use a rather small time step assuring 
l/A/3 ~ 2 * IQ^K. Also, the treatment of the interactions differs from our scheme: the 
authors of Ref. ||32| perform a numerical solution of the Bloch equation for the two-particle 



density matrix whereas we use an analytical approximation for the effective pair interaction 
(based on the Kelbg potential, see above). Finally, Ref. approximately computes the 
nodal surface of the density matrix using a variational ansatz which is then used to restrict 
the integrations to the region of positive density matrix. For more details regarding the 
restricted PIMC simulations, see Refs. [ p!9| , p^ . 
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FIG. 2. Comparison of direct and restricted PIMC results and analytical results (FADE) 
for the pressure and total energy of dense hydrogen as a function of temperature for = 1.86, 
corresponding to n = 2.5 • For illustration, also the coupling and degeneracy parameters 

F and ngA^ are shown in the upper figure. 

Let us now turn to a comparison of the numerical results. The restricted PIMC simulation 



data for dense hydrogen are taken from Ref. [32|. A comparison of results for the pressure 



and the internal energy for a fixed value of the density (r^ = 1.86) is shown in Fig. 1 and 
TABLE I. At high temperatures, above 50,000 K, where only a small fraction of atoms is 
expected, the agreement is rather good. This is remarkable since the nonideality and the 
degeneracy reach values of 3 and 10, respectively. This result demonstrates that, at least for 
~ 1 and for T > 50, OOOi^ both methods yield results which are more or less equivalent. At 
T < 50, OOOK, where partial ionization is expected, we still observe a reasonable agreement 
of both approaches, however, we see also that the differences start to grow. The reasons for 
that are manyfold. From our results we conclude that the main problem is not the bound 
state formation - atoms and molecules are well described by the two PIMC simulations 
which use a physical picture which does not involve any artificial distinction between free 
and bound electrons, e.g. [^. On the other hand, with growing degeneracy nA^, both 
PIMC methods become less reliable, and a detailed analysis, although being very desirable, 
will have to be based on more extensive calculations in the future. 
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Further we present in TAB. I also Fade results for the weakly nonideal region. We find 
good agreement with the FIMC results for T > 10^ K. Details on the method of these 
analytical calculations will be discussed in the next section. 



V. COMPARISON WITH ANALYTICAL APPROXIMATIONS FOR THE 
THERMODYNAMIC FUNCTIONS OF STRONGLY IONIZED DENSE PLASMAS 

In this section we give a comparison of the available data points from direct FIMC 
calculations with analytical estimates based on Fade approximations for strongly ionized 
plasmas PJ3^-p6|. The comparison concentrates on H-plasmas in a region in the density- 
temperature plane with the following borders 

0.2 < F < 1.6, 

0.2 < n^kl < 5, (21) 

which will be called "rhomb of moderate nonideality and moderate degeneracy" (see bold 
rhomb in Fig. 1). With respect to analj^ical treatment, this rhombic region is of particular 
difficulty since none of the known analytic limiting expressions is valid. Further we calculated 
several points for = 1.86 and F < 2 which correspond to the FIMC data discussed in 
the previous section and also an isotherm at T = 50, OOOK including some data at higher 
density, outside the rhomb, cf. Fig. 1 for an overviev. 

We demonstrate below that the Fade approximations which interpolate between the 
limits where theoretical results are available are a useful tool for the description of the 
available data points, at least for the case of moderate nonideality F < 1.6 and moderate 
degeneracy UeA^ < 5. The Fade approximations which we use here were constructed in 
earlier work, [Bl-|36|, from the known analytical results for limiting cases of low density 



||3|j33[| and high density 0. The structure of the Fade approximations was devised in such 
a way that they are analytically exact up to quadratic terms in the density (up to the 
second virial coefficient) and interpolate between the virial expansions and the high-density 
asymptotic expressions P^-|5B|. The formation of bound states was taken into account by 
using a chemical picture. This means the plasma is considered as a mixture of free electrons, 
free ions, atoms and molecules which are in chemical equilibrium, being described by mass 
action laws or minimization of the free energy . 

We follow in large here this cited work, only the contribution of the ion-ion interaction 
which is, in most cases, the largest one, was substantially improved following recent work of 
Kahlbaum, who succeeded in describing the available classical Monte Carlo data for the ions 
by accurate Fade approximations |^B|. By using Kahlbaum's formulas we achieve a rather 



accurate description of the thermodynamics in the region where the plasma behaves like a 
classical one-component ion plasma imbedded into a sea of nearly ideal electrons. This is 
the region where the electrons are strongly degenerate 

UeAl > 1 and < 1, (22) 

and the ions are still classical but nonideal 

F > 1 and UiA^ < 1. (23) 
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This region lies in the upper left corner of Fig. 1. 

With respect to the chemical picture we restrict ourselves to the region of strong ioniza- 
tion where the number of atoms is still relatively low and where the fraction of molecules 
is small as well, see below. We will discuss here only the general structure of the Fade 
formulae. For example, the internal energy density of the plasma is given by 



u 



Uid + Uint- (24) 



Here Uid is the internal energy of an ideal plasma consisting of Fermi electrons, classical 
protons and classical atoms, and Umt is the interaction energy 

Uint = Uii + Uee + Uie + UydW ■ (25) 

The interaction contribution to the internal energy consists of four terms: 

• Ion-ion interaction contribution: this term which, in general, yields the largest con- 
tribution is generated by the OCP subsystem of the protons. For the OCP energy of 
protons many expressions are available, e.g. |^5|. We have used here the most precise 
formula due to Kahlbaum |^B|] which interpolates between the Debye region, ua ~ F^/^, 



and the high density fluid, Uu ~ F. 

Electron-electron interaction: This term corresponds to the OCP energy of the electron 



subsystem. We used the rather simple expressions used in earlier work p^ , p5 



Electron-proton interaction: This term corresponds to the interaction between the two 
OCP subsystems which is mostly due to polarization effects. Again, we used the rather 
simple expressions proposed in earlier work p4| , p5| . 



• Van der Waals contribution: In the region of densities and temperatures defined above 
this contribution gives only a small correction. Therefore, this term was approximated 
here in the simplest way by a second virial contribution. The neutral particles were 
treated as hard spheres. 

In the region of densities which are studied here, molecules do not play a role, therefore, 
the formation of molecules was taken into account only in a very rough approximation 
according to Ref. The number density of the neutrals was calculated on the basis of a 



nonideal Saha equation. We restricted this comparison to a region where the number density 
of neutrals is relatively small, the degree of ionization being larger than 75%. 

The contributions to the pressure were calculated, in part, from scaling relations e.g. we 
used Pa = Uu/S, and, for the other (smaller) contributions, by numerical differentiation of 



the free energy given earlier ||3^j35[| . In a similar way, the chemical potential which appears 
in the nonideal Saha equation was obtained. For the partition function in the Saha equation 
we used the Brillouin-Planck-Larkin expression |]3|,|36|. The solution of the nonideal Saha 
equation which determines the degree of ionization (the density of the atoms) was solved by 
up to 100 iterations starting from the ideal Saha equation. 
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FIG. 3. Comparison of Fade calculations (lines without symbols) for the internal energy with 
the direct FIMC results (lines with full circles). 

Since all the expression described so far are given in analytic form, the calculation of 
about 1000 data points for energy and pressure takes less than a minute on a PC. The result 
of our calculations for density-temperature points in the "rhomb of moderate nonideality 
and moderate degeneracy" are given in Figs. 2,3. Further, we give in TAB. I several data 
points obtained from the Fade formulas. Since the Fade formulas used here do not apply to 
low temperatures, we included in TAB. I only Fade data for T > lO^i^. 
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FIG. 4. Comparison of Fade calculations (lines without symbols) of the pressure (in units of 
the Boltzmann pressure) with direct FIMC simulation results (lines with full circles). 

Summarizing the results for the internal energy and for the pressure we find that the Fade 
results, with a few exceptions, agree well with the FIMC data in the region of the density 
temperature plane, where F < 1.6 and nA^ < 5. The agreement is particularly good for the 
energies. [The larger deviations for the pressure may be due to the numerical differentiation.] 
In fact, the Fade formulas used here in combination with the chemical picture works only 
in the case that the plasma is strongly ionized, i.e. the degree of ionization is larger than 
75%. The description of the region where a higher percentage of atoms and, due to this, 
also molecules is present needs a more refined chemical picture [^^|8 . 
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FIG. 5. Comparison of Fade calculations (lines without symbols) of the pressure (in units of the 
Boltzmann pressure) with direct FIMC simulation results (lines with full circles) for an isotherm 
T = 50, OOOiC. 

Finally, we compare the Fade and FIMC data along the isotherm T = 50, OOOi^' which 
is given in Fig. 5. This figure shows the transition from a classical ideal gas (low density) 
to a nearly ideal quantum gas (limit of high density). In the central part, n < lO^^cm"^ < 
lO^^cm"^, Coulomb interaction leads to strong deviations from the behavior of an ideal 
plasma. The strong increase of the energy at high density is due to the Mott effect and to 
the increase of the ideal quantum contribution to the electron energy. Comparing the Fade 
and FIMC results, we find good agreement up to electron densities n = lO^^cm^. For higher 
densities, the deviations are growing. For intermediate densities, n < lO^^cm"^ < lO^'^cm"^ 
the PIMC data are more rehable. On the other hand, in the limit of very high density, -C 1, 
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the Pade results are known to correctly approach the ideal quantum plasma limit whereas the 
PIMC data should be regarded as preliminary due to the extremely high electron degeneracy. 
Interestingly, we find that at high density the Fade data approaches the ideal curves earlier 
than the PIMC data which is important for further improvement of the presented Monte 
Carlo approach. 

VI. DISCUSSION 

This work is devoted to the investigation of the thermodynamic properties of hot dense 
partially ionized plasmas in the pressure range between 0.1 and 100 Mbar. Most of the new 
results are based on a Quantum Monte Carlo study of a correlated proton-electron system 
with degenerate electrons and classical protons. In this paper, we gave a detailed derivation 
of improved estimators for the internal energy and the equation of state for use in direct 
fermionic path integral simulations. Also, we gave a rigorous justification for the use of an 
effective quantum pair potential (Kelbg potential) in PIMC simulations. 

Further, we compared our direct PIMC results with independent restricted PIMC data 
of Militzer and Ceperley for one isochor corresponding to = 1.86, Fig. 2. We found very 
good quantitative agreement between the two PIMC methods for temperatures in the range 
of 50, OOOK < T < lO^K, where F < 3 and ngA^ < 10. This region is particularly com- 
plicated as here pressure and temperature ionization occur and, therefore, an accurate and 
consistent treatment of scattering and bound states is crucial. This agreement is remarkable 
because the two simulation methods are completely independent and use essentially different 
approximations. We, therefore, expect that the results for the thermodynamic properties 
of high pressure hydrogen plasmas in this temperature-density range are reliable within the 
limits of the simulation accuracy. This is the main result of the present paper. 

In future work, it will be important to extend the range of agreement. To analyze 
the deviations between the simulation methods, we also included some data for = 1.86 
and lower temperatures, 10, OOOK < T < 50, 000-ft", Fig. 2. At this point, no conclusive 
answer about the reasons of the deviations can be given. For these parameters, the electron 
degeneracy is growing rapidly and, therefore, each of the simulation methods is becoming 
less reliable. So these data should be regarded as preliminary results which will be useful 
for future improvements of the simulations. 

Furthermore, the Monte Carlo results allowed us to develop and test analytical approxi- 
mations of Pade-type which are improvements of earlier approximations [P,P^-|5B| in a region 
in the density-temperature plane bounded by F < 1.6 and n^Al < 5. This is a region of 
moderate nonideality and degeneracy and high degree of ionization. We have shown that for 
these parameters, the Fade approximations which interpolate between the limits where the- 
oretical results are available agree well with the Monte Carlo data, cf. Figs. 2-4 and Table I. 
Thus, these approximations provide a useful tool for the description of these plasmas which 
include hydrogen at a pressur between 0.1 and 100 Mbar. At lower temperature, deviations 
from the Monte Carlo data are growing, cf. Fig. 2. This is mostly due to the growing 
role of bound states. Whether the Pade approximations, in combination with an improved 
chemical picture (mass action law), continue to work at lower temperatures, has still to be 
explored, first steps are under way [|^. 
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Also, we showed some data for T = 50, OOOK and higher pressure, up to p ~ 10^ Mbar, 
Fig. 5. Here the Monte Carlo simulations are particularly difficult due to the high electron 
degeneracy, and they can beneffi from the Fade simulations, as the latter correctly reproduce 
the high-density limit, <^ 1. 
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TABLES 



TABLE I. Direct versus restricted PIMC [23] simulation results (upper and middle lines, 
respectively) and results of Pade calculations (numbers in the lowest lines) for the pressure p(Mbar) 
and energy £^(2NRy) for dense hydrogen (deuterium [23]) for = 1.86 



T, lOOOK 




r 


p, Mbar 


E, 2NRy 


1000 


0.10 


0.169 


67.74 ± 0.02 
66.86 ± 0.08 
67.38 


9.050 ± 0.005 
9.018 ± 0.015 
9.063 


500 


0.29 


0.339 


32.85 ± 0.03 
32.13 ± 0.05 
31.91 


4.169 ± 0.003 
4.114 ± 0.007 
4.162 


250 


0.83 


0.679 


15.37 ± 0.01 
14.91 ± 0.03 
14.40 


1.654 ± 0.005 
1.629 ± 0.007 
1.679 


125 


2.33 


1.350 


6.98 ± 0.01 
6.66 ± 0.02 
6.47 


0.412 ± 0.005 
0.404 ± 0.004 
0.471 


62.5 


6.58 


2.701 


3.07 ± 0.02 
2.99 ± 0.04 


-0.248 ± 0.005 
-0.140 ± 0.007 


31.25 


18.48 


5.376 


2.20 ± 0.01 

1.58 ± 0.07 


-2.377 ± 0.005 

-0.360 ± 0.010 
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